from pylab import linspace, sin, pi, cos, exp, double

def hilite(string, status, bold=False):
    """Just to Add some color and fun"""
    attr = []
    if status:
        # green
        attr.append('33')
    else:
        # red
        attr.append('31')
    if bold:

        attr.append('1')
    return '\x1b[%sm%s\x1b[0m' % (';'.join(attr), string)


## Defining a point in two dimension ##
class Point():
    def __init__(self, x, y):
        self.x = double(x)
        self.y = double(y)
    def __add__(self,p):
        t=Point(0,0)
        t.x=self.x+p.x
        t.y=self.y+p.y
        return t
    def __sub__(self,p):
        t=Point(0,0)
        t.x=self.x-p.x
        t.y=self.y-p.y
        return t
    def __mul_(self,p):
        t=Point(0,0)
        t.x=self.x*p.x
        t.y=self.y*p.y
        return t
    def __abs__(self):
        return (self.x**2 + self.y**2)**.5
## Norm of vector in two dimension ##
def Norm(a, b):
    return ((a.x - b.x)**2+(a.y - b.y)**2)**.5



## Two dimensional gradient function return a two dimensional point ##
##def TwoDGrad(kernel, a, b,h):
    

class interpolatorSPH():
    """Content """

    reqpoints = linspace(-pi,pi,13)
    nodes = linspace(-2*pi,2*pi,400)
    fval = sin(nodes)

##Function Definition follows##

 ##Cubic Spline Kernel##

    def CubicSplineKernel(self, x0, x, h, *n):
        """ This is a Cubic Spline Kernal calculator for one dimension.
        Input: x0->Point at which Kernel is to be centered.
               x-->Point at which valve is to be calculated.
               h-->Step size.
        Output: Type<float>   
        """
        x0 = float(x0)
        x = float(x)
        h = float(h)
        q = abs((x0-x)/h);
        fac = ( 2/(3*h))
        if q <= 1 and q >= 0:
            val=fac*(1-((3*q*q)/2)+((3*(q**3))/4))
        elif q > 1 and q <= 2:
            val = fac*(((2-q)**3)/4)
        else:
            val = 0
        return val


    def gradCubicSplineKernel(self, x0, x, h, *n):
        """ Gives Derivative of the One Dimensional Cubic Spline Kernel defined above.
        Input:x0->Point at which Kernel is to be centered.
              x-->Point at which derivative of Kernel is to be calculated.
              h-->Step size.
        Output: type<float>
        """
        x0 = float(x0)
        x = float(x)
        h = float(h)    
        r = abs(x0-x)
        q = (r)/h
        fac = (2/(3*h))
        if r == 0:
            dqbydx = Point(0,0)
        else:
            dqbydx = Point( -(x0-x).x/(h*(r)), -(x0-x).y/(h*(r)))


        if q >= 0 and q <= 1:
            grad = Point(fac * ((-3.0 * q) + ((9.0 * q * q)/4)) * dqbydx.x, fac*((-3 * q) + ((9 * q * q)/4)) * dqbydx.y)
            
        elif q > 1 and q <= 2:
            grad = Point(-fac * (3 * ((2-q) * (2-q))/4)* dqbydx.x, -fac * (3 * ((2-q) * (2-q))/4) * dqbydx.y)

        else:
            grad = Point(0,0)


        return grad

 ##Cubic Spline Kernel Ends##

 ##    Wn Kernel           ##

    def WnKernel(self, x0, x, h, n=4):
        """ This is a W n kernel calculator for one dimension. As given in 
        A one-parameter family of interpolating kernels for smoothed particle hydrodynamics studies'
        
        Input: x0->Point at which value is to be centered
               x-->Point at which valve is to be calculated
               h-->Step size
               n-->Parameter
        Output: type<float>   
        """
        x0 = float(x0)
        x = float(x)
        h = float(h)
        q = abs((x0-x)/h);
        fac = self.wnKernelFactor(h,n)
        
        if q == 0:
            val = fac*1
        elif q > 0 and q <= 2:
            val = fac*((sin((pi*q)/2))/((pi*q)/2))**n
        else:
            val = 0
        return val
      
    def wnKernelFactor(self, h, n=4):
        """Calculate the Normalisation factor of W n Kernel 
        Input: h-->Step size
               n-->Parameter
        Output: type<float>
        """
        coeffs = [.2645649, .1824975, -0.02426267, .003112410, -0.0002404560, .000008032609]
        Kn = 0
        power = 1
        for i in coeffs:
            Kn+= i*power
            power = power*n
    
       
        fac = Kn/h
        return fac

    def gradWnKernel(self, x0, x, h, n=4):
        """Calculate the gradient of W n Kernel 
        Input: x0->Point at which kernel is to be centered
               x-->Point at which grad valve is to be calculated
               h-->Step size
               n-->Parameter
        Output: type<float)
        """
        
        x0 = float(x0)
        x = float(x)
        h = float(h)    
        r = abs(x0-x)
        q = (r)/h
        fac = self.wnKernelFactor(h,n)
        if r == 0:
            dqbydx = 0
        else:
            dqbydx = -(x0-x)/(h*(r))
        tmp = (pi*q)/2
        if q == 0:
            grad = 0
        elif q > 0 and q <= 2:
            grad = dqbydx*(n*((sin(tmp)/(tmp))**(n-1)))*(tmp*(pi/2)*cos(tmp)-(pi/2)*sin(tmp))/(tmp**2)
        else:
            grad = 0
        return fac*grad
 ##    Wn Kernel Ends!   ##

 ##M6 Spline Kernel##

    def M6SplineKernel(self, x0, x, h, *n):
        """ This is a M6 Spline Kernal calculator for one dimension.
        Input: x0->Point at which Kernel is to be centered.
               x-->Point at which valve is to be calculated.
               h-->Step size.
        Output: Type<float>   
        """
        x0 = float(x0)
        x = float(x)
        h = float(h)
        q = abs((x0-x)/h);
        fac = (243.0/(2560.0*h))
        if q <= 2.0/3 and q >= 0:
            val = fac*(((2-q)**5)-6*(((4.0/3)-q)**5)+15*(((2.0/3)-q)**5))
        elif q > 2.0/3 and q <= 4.0/3:
            val = fac*(((2-q)**5)-6*(((4.0/3)-q)**5))
        elif q > 4.0/3 and q <= 2.0:
            val = fac*((2-q)**5)
        else:
            val = 0

        return val


    def gradM6SplineKernel(self, x0, x, h, *n):
        """ Gives Derivative of the One Dimensional M6 Spline Kernel defined above.
        Input:x0->Point at which Kernel is to be centered.
              x-->Point at which derivative of Kernel is to be calculated.
              h-->Step size.
        Output: type<float>
        """
        x0 = float(x0)
        x = float(x)
        h = float(h)    
        r = abs(x0-x)
        q = (r)/h
        fac = (243.0/(2560.0*h))
        if r == 0:
            dqbydx = 0
        else:
            dqbydx = -(x0-x)/(h*(r))

        if q <= 2.0/3 and q >= 0:
            grad = -fac*dqbydx*(5.0*((2.0-q)**4)-30*(((4.0/3)-q)**4)+75.0*(((2.0/3)-q)**4))
        elif q > 2.0/3 and q <= 4.0/3:
            grad = -fac*dqbydx*(5.0*((2.0-q)**4)-30.0*(((4.0/3)-q)**4))
        elif q > 4.0/3 and q <= 2:
            grad = -fac*dqbydx*(5.0*(2.0-q)**4)
        else:
            grad = 0

        return grad

 ##M6 Spline Kernel Ends##

 ##Gaussian Kernel##

    def GaussianKernel(self, x0, x, h, *n):
        """ This is a Gaussian Kernal calculator for one dimension.
        Input: x0->Point at which Kernel is to be centered.
               x-->Point at which valve is to be calculated.
               h-->Step size.
        Output: Type<float>   
        """
        x0 = float(x0)
        x = float(x)
        h = float(h)
        q = abs((x0-x)/h);
        fac = (1.0/((pi**.5)*h))
        val = fac*exp(-q*q)
        return val


    def gradGaussianKernel(self, x0, x, h, *n):
        """ Gives Derivative of the One Dimensional Gaussian Kernel defined above.
        Input:x0->Point at which Kernel is to be centered.
              x-->Point at which derivative of Kernel is to be calculated.
              h-->Step size.
        Output: type<float>
        """
        x0 = float(x0)
        x = float(x)
        h = float(h)    
        r = abs(x0-x)
        q = (r)/h
        fac = (1.0/((pi**.5)*h))

        if r == 0:
            dqbydx = 0
        else:
            dqbydx = -(x0-x)/(h*(r))

        grad = -2*q*fac*exp(-q*q)*dqbydx

        return grad

 ##Gaussian Kernel Ends##

 ##Cubic Spline Kernel##

    def W10Kernel(self, x0, x, h, *n):
        """ This is a W10 Kernel calculator for one dimension.
        Input: x0->Point at which Kernel is to be centered.
               x-->Point at which valve is to be calculated.
               h-->Step size.
        Output: Type<float>   
        """
        x0 = float(x0)
        x = float(x)
        h = float(h)
        q = abs((x0-x)/h);
        coeffs=[0.676758,-0.845947,0.422974,-0.105743,0.0132179,-0.000660896]
        fac = (1/h)
        if q>=0 and q<=2:
            val=0
            power=0
            for i in coeffs:
                val+=i*(q**(power))
                power+=2
        else:
            val=0
        return val*fac


    def gradW10Kernel(self, x0, x, h,*n):
        """ Gives Derivative of the One Dimensional W8 Kernel defined above.
        Input:x0->Point at which Kernel is to be centered.
              x-->Point at which derivative of Kernel is to be calculated.
              h-->Step size.
        Output: type<float>
        """
        x0 = float(x0)
        x = float(x)
        h = float(h)    
        r = abs(x0-x)
        q = (r)/h
        if r == 0:
            dqbydx = 0
        else:
            dqbydx = -(x0-x)/(h*(r))
        coeffs=[0.676758,-0.845947,0.422974,-0.105743,0.0132179,-0.000660896]
        fac = (1/h)
        if q>=0 and q<=2:
            val=0
            power=2
            for i in coeffs[1:]:
                val+=i*power*(q**(power-1))
                power+=2
        if q>=2:
            val=0
        return val*fac*dqbydx
        

 ##W10 Kernel Ends##




 ## One Dimensional Interpolation Function ##    

    def OneDInterpolation(self,req,nodes,funval,ratio,kernel,*n):
        """This function interpolates function having value (funval at src) at the points given by the list req.
        input: req   ->List of points at which to Interpolate.
               nodes ->List of points at which function value is known.
               funval->Function value at points given by list nodes.
               ratio ->h/dx.
               test  ->Kernel to use.
               It is assumed that difference between two consective node values is constant.
        Output: type<list>
        """
        result = []
        dx = abs(float(nodes[1])-float(nodes[0]))
        h = ratio*dx
        for val_dst in req:
            value = 0
            counter = 0
            val_dst = float(val_dst)
            for val_src in nodes:
                val_src = float(val_src)
                w = kernel(val_dst, val_src, float(h),*n)
                value+= w*float(funval[counter])*dx
                counter+=1
            result.append(value)
        return result

 ## InterPolation function Ends ##

 ## Test Kernel ##

    def test_Kernel(self):
        """Composed of two function testWn and testCubicSpline.
           Perform the following test on both the kernel with default values of arguments:
               -Integrate value of kernel and gradient and display the result.
                Ideally integration of kernel should give 1 and integration of gradient should give 0.
               -Calculate value of kernel at distance -+2h from centre and display the result.
                Ideally they should be very small, close to zero.
               -Calculate the gradient at centre and display the result. Ideally it should be 0.
           For changing argument check the definition.
           Range=linspace(-8,8,1000)
           """
        t=linspace(-1,1,21)
        def testGaussian(t):
            print('')
            print(hilite('Testing Gaussian Kernel....',False,True))
            x0 = 0
            h = 2*(t[1]-t[0])
           
            gradintegral = 0
            valintegral = 0
            for i in t:
                gradintegral+= (self.GaussianKernel(x0, i, h))*(t[1]-t[0])
                valintegral+= (self.GaussianKernel(x0, i, h))*(t[1]-t[0])
            print 'x0=',hilite(x0,True), 'h=',hilite(h,True)
            print('Integral of Values')
            print(hilite(valintegral,True))    
            print('Integral of gradient')
            print(hilite(gradintegral,True))
            print('Value at x0+2h positive no.')
            print(hilite(self.GaussianKernel(x0,x0+2*h,h),True))
            print('Value at x0-2h negative no.')
            print(hilite(self.GaussianKernel(x0,x0-2*h,h),True))
            print('Gradient at center')
            print(hilite(self.gradGaussianKernel(x0,x0,h),True))
            print(hilite('Testing of Gaussian Kernel Ends!',True,True))

        def testCubicSpline(t):
            print('')
            print(hilite('Testing Cubic Spline Kernel....',False,True))
            x0 = 0
            h = 2*(t[1]-t[0])
            
            gradintegral = 0
            valintegral = 0
            for i in t:
                gradintegral+= (self.gradCubicSplineKernel(x0, i, h))*(t[1]-t[0])
                valintegral+= (self.CubicSplineKernel(x0, i, h))*(t[1]-t[0])
            print 'x0=',hilite(x0,True), 'h=',hilite(h,True)
            print('Integral of Values')
            print(hilite(valintegral,True))    
            print('Integral of gradient')
            print(hilite(gradintegral,True))
            print('Value at x0+2h positive no.')
            print(hilite(self.CubicSplineKernel(x0,x0+2*h,h),True))
            print('Value at x0-2h negative no.')
            print(hilite(self.CubicSplineKernel(x0,x0-2*h,h),True))
            print('Gradient at center')
            print(hilite(self.gradCubicSplineKernel(x0,x0,h),True))
            print(hilite('Testing of Cubic Spline Kernel Ends!',True,True))

        def testM6Spline(t):
            print('')
            print(hilite('Testing M6 Spline Kernel....',False,True))
            x0 = 0
            h = 2*(t[1]-t[0])
            gradintegral = 0
            valintegral = 0
            for i in t:
                gradintegral+= (self.gradM6SplineKernel(x0, i, h))*(t[1]-t[0])
                valintegral+= (self.M6SplineKernel(x0, i, h))*(t[1]-t[0])
            print 'x0=',hilite(x0,True), 'h=',hilite(h,True)
            print('Integral of Values')
            print(hilite(valintegral,True))    
            print('Integral of gradient')
            print(hilite(gradintegral,True))
            print('Value at x0+2h positive no.')
            print(hilite(self.M6SplineKernel(x0,x0+2*h,h),True))
            print('Value at x0-2h negative no.')
            print(hilite(self.M6SplineKernel(x0,x0-2*h,h),True))
            print('Gradient at center')
            print(hilite(self.gradM6SplineKernel(x0,x0,h),True))
            print(hilite('Testing of M6 Spline Kernel Ends!',True))
           
        def testWn(t):
            print('')
            print(hilite('Testing Wn Spline Kernel....',False,True))
            x0 = 0
            h = 2*(t[1]-t[0])
            gradintegral = 0
            valintegral = 0
            for i in t:
                gradintegral+= (self.gradWnKernel(x0, i, h))*(t[1]-t[0])
                valintegral+= (self.WnKernel(x0, i, h))*(t[1]-t[0])
            print 'x0=',hilite(x0,True), 'h=',hilite(h,True)
            print('Integral of Values')
            print(hilite(valintegral,True))    
            print('Integral of gradient')
            print(hilite(gradintegral,True))
            print('Value at x0+2h positive no.')
            print(hilite(self.WnKernel(x0,x0+2*h,h),True))
            print('Value at x0-2h negative no.')
            print(hilite(self.WnKernel(x0,x0-2*h,h),True))
            print('Gradient at center')
            print(hilite(self.gradWnKernel(x0,x0,h),True))
            print(hilite('Testing of Wn Kernel Ends!',True,True))

        def testW10(t):
            print('')
            print(hilite('Testing W10 Kernel....',False,True))
            x0 = 0
            h = 2*(t[1]-t[0])
            gradintegral = 0
            valintegral = 0
            for i in t:
                gradintegral+= (self.gradW10Kernel(x0, i, h))*(t[1]-t[0])
                valintegral+= (self.W10Kernel(x0, i, h))*(t[1]-t[0])
            print 'x0=',hilite(x0,True), 'h=',hilite(h,True)
            print('Integral of Values')
            print(hilite(valintegral,True))    
            print('Integral of gradient')
            print(hilite(gradintegral,True))
            print('Value at x0+2h positive no.')
            print(hilite(self.W10Kernel(x0,x0+2*h,h),True))
            print('Value at x0-2h negative no.')
            print(hilite(self.W10Kernel(x0,x0-2*h,h),True))
            print('Gradient at center')
            print(hilite(self.gradW10Kernel(x0,x0,h),True))
            print(hilite('Testing of W10  Kernel Ends!',True,True))
        testWn(t)
        testCubicSpline(t)
        testM6Spline(t)
        testW10(t)
        testGaussian(t)















if __name__ == '__main__':
    main()


 
    
